#ifndef EQUATIONSOLVER_HPP
#define EQUATIONSOLVER_HPP
#include <iostream>
#include "Function.hpp"
#include <cmath> 
#include <stdexcept>

inline int sgn(double x) {
    return (x > 0) - (x < 0);
}

class EquationSolver {
protected:
    const Function &F; 
public:
    EquationSolver(const Function& F) : F(F) {}
    virtual double solve() = 0; 
};

class Bisection_Method : public EquationSolver {
private:
    double a, b;
    double eps, delta;
    int max_iter;
public:
    Bisection_Method(const Function &F, double a, double b, 
                     double eps = 1e-7, double delta = 1e-6, int max_iter = 100)
        : EquationSolver(F), a(a), b(b), eps(eps), delta(delta), max_iter(max_iter) {}
    
    virtual double solve() override {
        if (sgn(F(a)) == sgn(F(b))) {
            throw std::runtime_error("Function values at the endpoints must have different signs.");
        }

        double h, u, c, w;
        int i;
        h = b - a;
        u = F(a);
        
        for (i = 0; i < max_iter; i++) {
            h /= 2; 
            c = a + h;
            if (std::abs(h) < delta || i == max_iter - 1) {
                break;
            }
            w = F(c);
            if (std::abs(w) < eps) {
                return c; 
            } else if (sgn(w) == sgn(u)) {
                a = c; 
            } else {
                b = c; 
            }
        }
        return c; 
    }
};

class Newton_Method : public EquationSolver {
private:
    double x0;
    double eps;
    int max_iter;
public:
    Newton_Method(const Function &F, double x0,
                  double eps = 1e-7, int max_iter = 100)
        : EquationSolver(F), x0(x0), max_iter(max_iter), eps(eps) {}
    
    virtual double solve() override {
        double x = x0;
        double fx = F(x);
        double dfx = F.derivative(x);
        for (int i = 0; i < max_iter; i++) {
            fx = F(x);
            dfx = F.derivative(x);
            if (std::abs(fx) < eps) {
                break; 
            }

            if (dfx == 0) {
                std::cout << "Derivative is zero; no solution found." << std::endl;
                break;
            }

            x = x-fx / dfx; 
        }
		return x;
    }
};

class Secant_Method : public EquationSolver {
private:
    double x0,x1;
    double eps,delta;
    int max_iter;
public:
    Secant_Method(const Function &F, double x0, double x1,
                  double eps = 1e-7,double delta = 1e-6,int max_iter = 100)
        : EquationSolver(F),x0(x0),x1(x1), eps(eps), delta(delta), max_iter(max_iter){}
    
    virtual double solve() override {
        double xx = x1;
        double x = x0;
        double u = F(xx);
        double v = F(x);

        for (int k = 2; k <= max_iter; ++k) {
            double s = (xx - x) / (u - v);
            x = xx;
            xx = xx - u * s;
            if (std::abs(xx - x) < delta) {
                return x; 
            }
            v = u;
            u = F(xx); 
            if (std::abs(u) < eps) {
                return xx;
            }
        }
        throw std::runtime_error("Max iterations reached.");
        return xx;
    }
};
#endif 

